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Least  Squares  Computations  in  Science  and  Engineering 

Abstract 

Least  squares  computations  constitute  a  fundamental  tool  in  science  and  engineering.  The  reason  is 
diat  they  play  a  critical  role  in  fitting  numerical  models  to  real  world  observations.  This  AFOSR 
supported  research  effort  has  been  concerned  with  the  design  and  testing  of  new  algorithms  for 
least  squares  computations  and  optimization  in  science  and  engineering.  The  objectives  were  to 
mathematically  develop,  test,  and  analyze  fast  numerical  algorithms  for  the  efficient  soluticm  to 
problems  on  modem  high  performance  computers.  The  focus  of  this  project  was  the  tq^plication  of 
scientific  computing  technology  in  the  area  of  signal  and  image  processing.  Very  many  problems 
lead  to  over  determined  systems  of  linear  or  nonlinear  equations  that  are  often  solved  by  least 
squares  or  related  t^timization  methods.  Generally,  tte  problems  are  accmnpanied  by  constraints, 
such  as  bound  constraints,  and  the  observations  are  comipted  by  noise.  The  project  has  involved 
the  application  of  scientific  computing  in  the  area  of  computational  linear  and  nonlinear  least 
squares  methods  with  particular  applications  in  image  and  signal  processing,  where  recovering 
images  is  often  an  ill-posed  inverse  problem.  Additional  work  included  control  computations 
associated  witii  adaptive  cities.  General  topic  areas  that  were  considered  in  this  research  include: 

•  Fast  algorithms  for  adaptive  estimation  and filtering,  along  with  sensitivity 
analysis,  with  2q)plications  in  signal  processing. 

Image  reconstruction  algorithms  for  3-D  computed  tomography,  in  the  veiificatirm 
and  inspection  of  aerospace  structures. 

•  Fast  deconvolution  methods  using  the  FFT,  including  regularization  techniques,  for 
certain  structured  ill-posed  inverse  problems  arising  in  2-D  image  restoration. 

•  yiasm  trace  maximization  appUcations  for  control  with  applications  to 

adaptive  optics  systems. 

Some  of  diese  research  projects  have  been  investigated  in  collaboration  witii  researchers  at  Phillips 
Air  Force  Laboratory,  Kirkland  AFB. 

Key  Words:  constrained  least  squares,  adaptive  filtering,  adaptive  optics,  deconvolution,  image 
resmration,  parallel  algorithms,  trace  maximization,  inverse  problems,  FFT. 
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1.  Introduction 

This  AFOSR  supported  research  effort  has  been  concerned  with  the  design  and  testing  of  new 
algorithms  for  least  squares  computations  and  optimization  in  science  and  engineering.  The 
objectives  were  to  mathematically  develop,  test,  and  analyze  fast  numerical  algorithms  for  the 
efficient  solution  to  scientific  problems  on  modem  high  performance  computers.  The  focus  of  dus 
project  was  application  of  scientific  computing  technology  in  the  area  of  signal  and  image 
processing. 

Special  emphasis  is  given  to  applications  of  potential  interest  to  die  U.S.  Air  Force.  In  particular, 
die  principle  investigator  visited  and  interacted  with  researchers  in  the  Lasers  and  Imaging  Division 
at  die  Phillips  Air  Force  Laboratory  at  Kirkland  AFB  on  three  occasions  (See  Section  3).  Fte  also 
^nt  the  first  six  months  of  1992  visiting  the  Institute  for  Mathematics  and  Its  Applications  at  the 
University  of  Minnesota.  The  visit  enabled  him  to  interact  and  collaborate  with  several  researchers 
fiom  universities,  industry,  and  government  laborannies  fiom  around  the  country.  During  the 
three  years  of  the  grant,  two  Ph.D.  students  were  directed  and  one  is  in  progress  (See  Section  4), 
twenty-four  papers  have  been  completed  or  are  currently  in  progress  (See  Section  5),  and  twenty- 
eight  technical  presentations  have  been  given  on  diis  AFOSR  ^nsmed  research  (See  Section  6). 
We  turn  now  to  a  brief  summary  of  the  research  objectives. 

Part  of  ourresearch  effort  was  in  the  area  time-recursive  least  squares.  One  of  the  main  efforts 
here  has  been  directed  toward  determining  stability,  conditioning,  perfcmnance,  and  architectural 
inqilementation  of  fast  adaptive  algorithms  for  FFT-based  time  recursive  methods.  Such  fast  algo¬ 
rithms  are  receiving  ever  increasing  attention  due  to  die  expanding  availability  (rf  commercial  multi¬ 
processors  and  special  purpose  computers.  The  results  of  this  effort  is  widely  ^plicable  to  the 
next  generatitm  of  signal  processing  systems  in  business,  defense,  and  engineering.  This  effnt  has 
invdved  interaction  with  researchers  in  closed-loop  active  noise  (vibration)  control  at  Phillips  Air 
FOTce  LabOTatory,  Kirkland  AFB. 

In  the  area  of  large  scale  least  squares  computations,  with  applications,  e.g.,  to  image 
reconstruction  and  restoration^  we  have  been  concerned  with  the  development  of  effective  FFT- 
based  preconditioners  for  linear  and  non-linear  conjugate  gradient  iterations.  Special  attention  was 
given  to  deconvolution  problems,  which  lead  to  block  Toeplitz  and  related  least  squares  problems. 
These  problems  arise  in  a  variety  of  application  areas,  especially  reconstruction  and  restoratitm 
medrads  in  image  processing.  An  important  problem  that  has  been  addressed  here  is  diat  of  image 
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restorati(Mi  of  the  object  from  a  degraded  and  noisy  inoage  measurement  This  effort  has  also 
involved  interaction  with  researchers  at  Phillips  Air  Force  Laboratory  at  Kirkland  AFB. 

Another  direction  of  our  work  involves  the  design  and  testing  of  scalable  parallel  algmithms  for 
signal  and  image  processing  on  massively  parallel  OHnputers.  This  work  is  in  conjunction  with 
personnel  at  General  Electric  CRD  in  Schnectady,  NY  and  at  Cray  Research  in  Minnellis,  MN. 
The  principle  investigatcn'  is  working  on  these  {noblems  with  my  fomer  Ph  J>.  student  William 
Harrod,  at  Cray  Research.  Technology  transfer  to  industry  for  our  algorithms  and  software  is  an 
important  aspect  of  this  work. 

Adaptive  qptics  control  systems  form  the  subject  of  considerable  recent  investigation.  This  project 
relates  to  a  very  interesting  problem  posed  Dr.  Brent  EUerbroek  during  one  of  the  principle 
investigat(n*s  visits  to  Hiillips  Laboratxvy,  Kirklaiui  AFB.  Specially  designed  deformable  minors 
(q>erating  in  a  closed-loop  adaptive  optics  system  can  partially  compensate  for  the  effects  of 
atmospheric  turbulence.  The  systems  detect  the  atnoospheric  distortions  using  either  a  natural  guide 
star  (point)  image  or  a  guide  star  artificially  generated  fiom  the  back  scatter  of  a  laser  generated 
beacon.  A  wave  front  sensor  measures  the  q)tical  distortions  that  can  then  be  partially  nullified  by 
deforming  a  flexible  minor  in  the  telescope.  The  problem  reduces  to  trace  maximization 
cooqtutations.  The  optimal  control  system  may  be  used  in  conjunction  with  the  Department  of 
Defnise's  largest  optical  telescope,  currently  being  installed  at  Kirkland  Air  Force  Base's  Starfire 
Optical  Range. 

General  u^c  areas  investigated  in  tiiis  research  iiKluded: 

•  Fast  algorithms  for  adaptive  estimadon  and  filtering,  along  with  sensitivity  analysis,  witii 
i^rplications  in  signal  processing. 

•  hnage  reconstnicticxi  algorithms  fix  3-D  computed  tomography,  in  die  verification  and 
inspectkm  of  aerospace  structures. 

•  Fast  deconvolution  methods  using  the  FFT,  including  regularization  techniques,  for  certain 
structured  ill-posed  inverse  problems  arising  in  2-D  image  restoration. 

•  Masnx  trace  maximization  (plications  for  control  /nerAodr,  with  applications  to  ad^ve 
qrtics  ^sterns. 

Specific  research  accoirq)lishments  and  activities  are  outlined  in  the  following  sections. 
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2,  Research  Accomplishments 


Experience  in  mathematical  analysis  and  scientific  computing  is  being  used  to  attack  some  difficult 
cmnputational  problems,  arising  especially  in  multidimensional  signal  and  image  processing.  The 
principle  investigator  has  work  on  the  implied  aspects  of  this  research  with  smne  colleagues  and 
former  Ph.D.  students  in  industrial  laboratories,  including,  Boeing,  Cray  Research,  and  General 
Electric,  as  well  as  with  some  researchers  widi  the  Air  Fence  Phillips  Laborattny,  Kirkland  AFB. 

♦  FFT-Based  Recursive  Least  Squares  Computations  with  Applications  to 
Adaptive  Filtering.  In  1986,  Gilbert  Strang  addressed  the  question  of  whether  iterative 
m^ods  can  compete  with  direct  methods  for  solving  symmetric  positive  definite  Toeplitz  systems 
of  linear  equations.  The  answer  has  turned  out  m  be  an  unqualified  yes.  Strang  proposed  the  use 
of  circulant  matrices  to  precondition  conjugate  gradient  iterations  for  Toeplitz  systems.  The  reason 
this  iqrproach  is  coixqwtitive  with  direct  methods  is  clear.  The  use  of  circulant  prectmditioners  for 
these  problems  allows  the  use  of  Fourier  transforms  throughout  the  computations,  and  these  FFT- 
based  iterations  are  not  only  numerically  efficient,  but  also  hi^y  parallelizaUe. 

Numerous  articles  have  extended  the  Strang  idea  to  more  general  Toeplitz  systems,  and  several 
types  of  circulant  preconditioners  have  been  suggested.  Recently,  we  have  developed  FFT-based 
preconditioners  for  solving  Toeplitz  least  squares  problems,  where  applications  include  such 
important  signal  processing  problems  as  active  noise  cancellatitm,  seismic  deconvolution,  data 
ctwapression,  and  image  restoration.  The  purpose  of  part  of  this  work  is  to  address  the  question  of 
whether  iterative  methods  can  compete  with  direct  methods  for  recursively  solving  least  squares 
proUems  in  an  adiq)tive  signal  processing  envirexunent. 

Adaptive  finite  impulse  response  (FIR)  filters  ate  used  extensively  in  many  signal  processing  and 
ccMittol  applications:  for  instance,  in  system  identification,  equalization  of  telephone  channels, 
spectrum  analysis,  noise  cancellation,  echo  cancellation  and  in  linear  predictive  coding.  The  main 
concerns  in  the  design  of  ad^tive  filtering  algorithms  are  tiieir  convergence  performance  and  their 
computaticHial  requirements.  These  concerns  are  especially  important  when  the  filters  are  used  in 
real-time  signal  processing  iq)plications  at  the  sizes  of  the  filters  are  very  large  (as  is  the  case  in 
acoustic  echo  or  active  noise  cancellation  problems.)  We  are  considering  both  recursive  least 
squares  (RLS)  computations  and  least  mean  squares  (LMS)  methods.  For  the  LMS  case,  we  have 
proposed  a  new  FFT-based  LSM-Newton  algorithtiL  Preliminary  numerical  results  show  that  the 
performance  of  the  algorithm  is  good  in  the  sense  of  stability  and  convergence.  Also,  the 
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complexity  of  the  algorithm  fw  n  filter  parameters  is  only  0(n  log  n)  operations  per  each  adaptive 
time  step  and  the  algorithm  itself  is  highly  parallelizable.  These  attractive  features  could  lead  to  the 
use  of  the  algorithm  in  diverse  ads^tive  filtering  iq)plications.  (See  publication  21.) 

Recursive  least  squares  problems  arise  naturally  in  uitoregressive  (AR)  and  AR  moving  average 
(ARMA)  computations  in  systems  identificadon,  and  in  general  FIR  filtering.  Thus  far  we  have 
reduced  the  classical  complexity  number  for  recursive  least  squares  computations  firom  O(n^)  to 
0(n  log  n),  by  using  the  FFT,  where  n  is  the  number  of  filter  parameters  (See  publicatimis  2, 
5,  6,  7,  8,  12,  13,  16,  19  and  21).  We  are  now  attempting  to  reduce  the  complexity  of  our 
transform-based  schemes  to  0(n),  by  incorporating  wavelet  iq>proximations.  Possible  applications 
iiKlude  adaptive  noise  cancellation  or  vibration  control  work  at  Phillips  Laboratory. 

♦  Scalable  Parallel  Algorithms  for  Image  Analysis  and  Reconstruction.  This  is 
part  of  some  3-way  collaboration  between  myself,  working  under  this  AFOSR  grant,  and 
personnel  at  the  General  Electric  CRD  at  Schnectady,  NY,  together  with  a  former  PhJ).  student, 
William  Harrod,  at  Cray  Research  in  Minneapolis.  Dr.  Harrod  is  in  charge  of  designing 
applications  software  for  Cray's  new  massively  parallel  computer,  their  Tourus  connection  3- 
dimensionai  architecture  (T3D). 

This  coUabmative  project  performed  by  the  private  industry  and  university  team  is  establishing  a 
scalable  parallel  programming  environment  for  message-passing  distributed  memray  architectures 
for  applications  in  multidimensional  signal  processing  and  image  analysis.  The  regularities  of  the 
data  structures,  the  homogeneous  nature  of  the  calculations,  and  the  number  of  applications  of 
multidimensional  signal  processing  and  image  analysis  provide  an  unparalleled  opportunity  to 
exploit  the  potential  of  modem  massively  parallel  computers.  We  are  concentrating  cm  develt^ing 
portable,  scalable,  parallel  algorithms  for  processing  and  analyzing  very  large  multidimensional 
data  sets  with  well-defined  neighborhood  relatitm^ps  between  data  elements.  Examples  of  such 
data  sets  include  X-ray  digital  radiogrtq>hs.  X-ray  computerized  tonoographs,  synthetic  aperture 
radar,  visual/uifirared  camera,  laser  radar,  towed  arrays,  and  high-resolution  satellite  images.  The 
algorithms  will  be  implemented  on  a  selection  of  target  architectures,  included  mesh  connected, 
hypercube,  heterogeneous,  and  re  configurable  arrays.  One  such  platform  includes  the  Cray  MPP 
siq)er  ctm^rater  T3D,  made  available  by  the  principle  investigaun's  former  Ph.D.  student  William 
Harrod  at  Cray  Research. 

Another  ttqnc  in  this  overall  project  is  the  implication  of  parallel  algorithms  to  image  reconstruction 
on  massively  parallel  architectures.  The  3-D  image  reconstraction  approach  we  are  considering 
with  General  Electric  is  based  on  the  Allowing: 
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•  Data  collection  by  cone  beam  scanning.  The  use  of  fan  or  parallel  beam  scanning  is  too  slow 
for  the  riata  collection  in  the  applications  we  have  in  mind. 

•  Conversion  of  the  cone  beam  projection  data  to  Radon  data,  consisting  of  planar  integrals. 

•  3-D  Radon  transform  inversion.  This  particular  method  consists  of  two  steps,  each  of  which 
involves  a  2-dimensional  filtered  back  projectitxi  algorithm. 

Some  feel  that  cone  beam  x-ray  computerized  tomogny>hy  (CT)  inspection  systems  will  form  the 
basis  for  the  next  generation  of  CT  inspection  technology  for  both  industrial  and  military  uses. 
However,  die  cone  beam  approach  is  still  too  computationally  intensive  for  many  applications  even 
on  advanced  parallel  architectures,  since  for  practical  problems  rapid  inspection  of  large  numbers 
of  objects  with  high  voxel  dimensions  is  needed;  tiius  the  interest  in  investigating  the  possible  use 
of  scalable  parallel  algorithms  for  MPP  architectures  hoe,  in  order  to  reduce  the  on-line  inspection 
titrte.  Contacts  through  General  Electric  are  have  also  been  made  with  researchers  at  the  Wright- 
Patterson  Laboratory,  concerning  our  help  with  their  project  of  periodic  non-destructive  inspection 
of  Air  Force  equipment  parts  and  materials.  (See  publications  10, 17  and  20.) 

♦  Deconvolution  Methods  for  Image  Restoration.  Here,  the  work  involves 
tqrplication  of  our  new  FFT-based  preconditioning  schemes  for  Toeplitz  and  related  least  squares 
computations  to  image  deblurring  problems.  Our  novel  approach  is  to  precondition  in  the 
frequency  domain,  while  iterating  in  the  spatial  domain  for  high  performance.  It  is  important  to 
recover  the  information  in  a  blurred  image;  in  remote  sensing,  details  about  the  photographed 
terrain  should  be  clarified;  in  medical  imaging,  the  diagnosis  is  based  on  the  clarity  of  the  x-ray 
radiographs  taken;  in  space  activities,  images  transmitted  to  earth  by  unmanned  or  manned 
spacecraft  need  to  be  analyzed.  Our  woik  on  prectmditioning  techniques  can  be  applied  directly  to 
image  restoration  computations  for  2-D  signals,  where  the  point  spread  function  (PSF), 
represented  by  h(Xyy),  is  spatially  invariant,  and  thus  leads  to  a  block  Toeplitz  operator.  Such 
situations  arise  in  ground  based  astronomical  imaging  tests  at  the  Startire  Optical  Range,  Kirkland 
AFB.  Equations  relating  the  observed  and  blurred  image  to  the  true  image  can  be  written 
algebraically  in  matrix  equation  form  as 


g  =  Hf +  T1, 

where  g  represents  the  blurred  image,  H  represents  the  discretized  PSF,  f  is  the  true  image,  and 
11  represents  the  noise  process  in  recording  the  image.  Practical  problems  are  typically  ill-posed 
and  therefore  very  sensitive  to  model  uncertainty  and  to  noise.  We  are  interacting  with  personnel  in 
the  Lasers  and  Imaging  Division  at  Phillips  Air  Force  Laboratory  on  this  project  (See  Section  3 
and  publications  10, 14, 17, 18,  and  24.) 
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♦  An  Optimization  Problem  in  Adaptive  Optics  Control.  This  work  was  not  part  of  the 
original  proposal,  but  came  about  when  an  interesting  research  problem  was  posed  during  a  visit 
by  the  principle  investigator  to  the  Hullips  Laboratory  in  February  1993,  by  Dr.  Brent  EUerbroek. 
The  <^>timization  {plication  involves  the  solution  of  a  trace  maximizatitm  |»oblem,  associated  with 
orthogonal  transformations  of  a  set  of  mirror  matrices  for  closed  loop  system  in  adq)tive  optics. 
This  adaptive  optics  work  relates  to  the  adaptive  control  of  a  set  of  very  fast-acting  deformable 
mirrors  whose  surface  shapes  can  change  rapidly  to  cmrect  for  optical  distortions  caused  by  the 
atmosphere.  Matrix  computations  often  play  an  important  role  in  aero-optics  applications  for 
astrtmomical  and  aero-optics  imaging.  The  problem  we  have  considered  involves  the  optimal  real¬ 
time  control  of  a  very  fast-acting  deformable  mirror  designed  for  atmospheric  turbulence 
compensation.  The  surface  shape  of  this  mirrors  must  change  rapidly  to  correct  for  time-varying 
optical  distortions  caused  by  the  atmosphere.  The  sensor  measurements  used  to  compute  these 
corrections  are  noisy,  and  for  each  spatial  mode  of  the  distortion  there  is  an  optinnal  control 
bandwidth  balancing  residual  errors  induced  by  sensOT  noise  and  servo  lag.  One  formulation  of 
this  problem  yields  a  functional  f(U)  to  be  maximized  over  unitary  matrices  U,  where  U  defines 
the  basis  of  orthonormal  control  noodes  to  be  used,  and  characterizes  adaptive  optics  performance 
when  aU  deformable  mirror  degrees  of  freedom  are  controlled  at  a  common  bandwidth.  Numerical 
ejqperiments  are  undertaken,  using  some  practical  data  fiom  an  astroncxnical  telescope  system  at  the 
Kirkland  AFB,  Starfire  Optical  Range.  The  problem  is  being  approached  using  eigenanalysis  (See 
Section  3,  and  publicatitms  22  and  23). 
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3.  Interaction  with  Phillips  AF  Laboratory 


The  focus  of  this  interaction  with  Phillips  Laboratory  at  Kirkland  AFB,  which  was  supported  by 
the  grant  AFOSR-91-0163,  is  the  ^^)plicadon  of  scientific  computing  in  the  area  of  computational 
methods  and  analysis  in  image  and  signal  processing.  General  topic  areas  that  have  been 
ctxisidered  in  afreet  of  the  AFOSR  research  project  include; 

1 .  Matrix  trace  nuvdmization  implications  to  control  of  defr^mable  mirrors  for  closed  loop 
adaptive  optics  systems,  in  conjunction  with  Brent  Ellerbroek  at  Phillips  Laboratory. 

2.  Fast  deconvolution  methods  using  FFT  and  fast  wavelet  based  conjugate  gradient 
iterations,  including  regularization  techrtiques  foe  ill-posed  inverse  problems  in  image 
resten^on,  in  conjunction  with  Richard  Gureras  and  Brent  Ellerbroek  at  Phillips. 

These  projects  are  briefly  reviewed  below.  Each  topic  relates  to  imaging  through  atmospheric 
turbulence. 

Astronomers  and  other  scientists  have  long  sought  to  overcome  the  degradation  of  astronomical 
image  quality  caused  by  the  effects  of  atmospheric  turbulence.  These  effects  are  in  part  due  to  the 
mixing  of  warm  and  cold  air  layers.  The  resulting  twinkling  of  the  stars  and  other  effects  are  the 
main  limitations  of  ground-based  imaging  for  both  scientific  and  defense  purposes.  As  pointed  out 
vividly  in  a  January  1994  National  Geographic  Magazine  article.  New  Eyes  on  the  Universe, 
exciting  technological  breakthroughs  are  rapidly  coming  to  the  aid  of  scientists  attempting  to  deblur 
astroncnnical  images. 

The  improvement  in  ground-based  image  quality  is  now  generally  attempted  in  two  stages.  The 
first  stage  occurs  as  the  observed  image  is  initially  formed.  Specially  designed  deformable  mirrors 
operating  in  a  closed-loop  adaptive  optics  system  can  partially  compensate  for  the  effects  of 
atmospheric  turbulence.  The  systems  detect  the  distorticxis  using  either  a  natural  guide  star  (point) 
image  or  a  guide  star  artificially  generated  from  the  back  scatter  of  a  laser  generated  beacon.  A 
wave  fiont  sensor  measures  the  optical  distortions  that  can  tiien  be  partially  nullified  by  deforming 
a  flexible  mirror  in  the  telescope.  To  be  effective,  these  corrections  have  to  be  performed  at  near 
real-time  speed.  Adaptive  optics  control  systems  form  the  subject  of  considerable  recent 
investigaticMi.  Project  #  1  relates  to  a  very  interesting  problem  posed  by  Dr.  Ellerbroek  during  one 
of  the  principle  investigator’s  visits  to  Phillips  Laboratory.  The  problem  relates  to  the  adaptive 
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control  of  a  set  of  very  fast-acting  deformable  mirrors  whose  surface  shapes  can  change  rapidly  to 
correct  for  optical  distortions  caused  by  the  atmosphere.  By  developing  a  mathematical  model  of 
the  closed  loop  adaptive  optics  system,  the  problem  of  estimating  the  optimal  mirror  control 
parameter  vectcn*  c(t)  at  time  t  has  been  reduced  by  Dr.  Ellerbroek  to  the  solution  of  a  matrix 
trace  maximization  problem,  associated  with  orthogonal  transformations  of  set  of  mirror  matrices 

Ml,  M2 . Mp,  for  the  closed  loop  system.  In  conjunction  with  C.  Van  Loan  and  N.  Pitsianis 

at  Cornell  University,  the  problem  has  been  approached  using  eigenanalysis.  Two  joint  papers 
with  Dr.  Ellerbroek  will  follow  finom  this  study.  The  optimal  control  system  may  be  used  in 
conjunction  with  the  Department  of  Defense's  largest  optical  telescope,  currently  being  installed  at 
Kirkland  Air  Force  Base's  Starfire  Optical  Range.  (See  publications  22  and  23.) 


The  second  stage  of  compensating  for  the  effects  of  atmo^heric  turbulence  generally  occurs  off¬ 
line,  and  consists  of  the  image  processing  step  of  restoration.  An  image  partially  corrected  by  the 
adaptive  optics  procedure  discussed  above  can  generally  be  enhanced  further  by  off-line  computer 
image  restoration.  Here,  large-scale  computations,  again  using  either  a  natural  guide  star  (point) 
image  or  a  guide  star  artificially  generated  from  the  back  scatter  of  a  laser  generated  beacon,  are 
used  to  deconvolve  the  blurring  effects  of  atmosph^c  turbulence.  In  this  regard,  removing  a 
linear,  shift-invariant  blur  from  a  signal  or  image  can  be  accomplished  by  inverse  or  Wiener 
Altering,  or  by  an  iterative  least  squares  debluning  procedure.  Because  of  the  ill-posed 
characteristics  of  the  deconvolution  problem,  in  the  presence  of  noise,  direct  filtering  methods 
often  yield  poor  results;  but,  on  the  other  hand,  iterative  methods  often  suffer  from  slow 
convergence  at  high  spatial  frequencies.  To  accelerate  convergence,  we  use  the  preconditioned 
conjugate  gradient  iterative  algorithm.  A  new  approximate  inverse  Toeplitz  preconditioner  is  used 
to  further  increase  the  rate  of  convergence.  Test  results  are  reported  for  a  ground-based 
astronomical  imaging  problem  in  publication  24.  In  particular.  Project  #2  is  concerned  with 

Here,  the  work  involves  application  of  our  new 
FFT-based  preconditioning  schemes  for  constrained  nonlinear  least  squares  computations  to  image 
debluning  problems.  Our  novel  approach  is  to  precondition  in  the  frequency  domain,  while 
iterating  with  nonlinear  conjugate  gradients  in  the  spatial  domain  for  high  performance.  This 
facilitates  application  of  constraints  such  as  non  negativity  and  rinite  support  Fast  throughput  is 
achieved  by  using  die  FFT  throughout  the  computations.  This  work  on  preconditioning  techniques 
can  be  applied  directly  to  image  restoration  computations  for  2-D  signals,  where  the  point  spread 
function  (PSF)  is  ^tially  invariant  Here  the  blurred  image  fnmation  process  is  modeled  as: 


g(x,y)  =  /  J  h(x-o,y-p)  f(a,p)  dodp  +  Tl(x,y). 


•oe  *00 
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Equations  relating  the  bluned  image  to  the  true  image  can  thus  be  written  in  discrete  matrix  form  as 
g  3  Hf  +  Ti,  where  g  represents  the  blurred  image,  H  represents  the  discretized  PSF,  f  is  the 
true  image,  and  T)  represents  the  noise  process  in  recording  the  image.  In  actual  practice,  no 
matrices  need  to  be  formed  to  execute  the  algorithm.  Image  degradation  occm  in  a  variety  of 
imaging  systems;  including  industrial  radiography,  non-destmctive  evaluation,  medical  imaging, 
remote  sensing,  aerial  reconnaissance,  ground-based  imaging,  and  space-based  imaging.  Practical 
problems  are  typically  ill-posed  and  therefore  very  sensitive  to  model  uncertainty  and  to  noise.  Our 
computations  are  being  performed  on  a  vector-parallel  Cray  Y-MP.  For  the  important  case  where 
the  point  spread  function  is  not  spatially  invariant  due  to  the  small  size  of  the  isoplanatic  angle  in 
atmospheric  imaging  systems,  the  use  of  wavelet-based  preconditioners  is  also  being  investigated. 
(See  publicaticm  20.) 

Interactitm  on  this  image  restoration  work  has  been  made  with  researchers  at  Phillips  LabOTatory  at 
Kirkland  AFB.  In  particular,  we  are  testing  our  schemes  on  data  used  at  Kirkland  from  an 
ensemble  of  photographs  taken  by  the  Air  Force  of  satellites  from  the  ground  looking  up.  The 
blurring  is  caused  by  atmospheric  ttirbulence.  This  area  of  research  may  lead  to  important 
advances  in  high  speed  ctxnputations  in  restoration  of  ground-based  images.  We  are  coordinating 
our  work  with  Brent  EUerbroek  and  with  Richard  Carreras,  who  are  kindly  providing  us  with  the 
real  imaging  data  at  Phillips  Laboratory  on  this  project  In  particular,  the  figure  below  illustrates  the 
fast  convergence  of  our  FFT-based  preconditioned  conjugate  gradient  method  for  data  provided  to 
us  from  Phillips  Labwaioty.  (The  figure  will  appear  in  publicatitm  24.) 
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4.  Graduate  Students _ 

The  following  Ph.D.  students  have  worked  with  or  are  currently  working  with  the  principal 
investigator  as  part  of  this  research  grant.  Two  of  the  students  graduated  at  the  end  of  1991,  and  a 
third  student  has  begun  her  Ph.D.  studies  in  G>mputer  Science  at  North  Carolina  State  University 
under  the  direction  of  the  principal  investigator. 

•  James  Nagy.  Dr.  Nagy  completed  his  Ph.D.  in  Mathematics  with  a  minor  in  Electrical  and 
Computer  Engineering  in  December  1991,  under  the  direction  of  the  principal  investigator.  He  was 
supported  for  (me  summer  under  this  grant  His  dissertation  topic  at  NC  State  University  involved 
Toeplitz  and  related  least  squares  computations  in  signal  and  image  processing,  identification,  es¬ 
timation  and  control.  He  was  awarded  a  Post  Doctoral  position  with  the  Institute  for  Mathematics 
and  Applications,  University  of  Minnesota,  for  the  academic  year  1991-92.  He  accepted  a  faculty 
posititm  in  Mathematics  at  Southern  Methtxlist  University  beginning  August  1992. 

•  William  Ferng.  Dr.  Femg  completed  his  Ph.D.  in  Mathematics  with  a  minor  in  Computer 
Science  in  December  1991,  under  the  direction  of  the  principal  investigator.  He  was  supported  for 
one  summer  under  this  grant  He  did  his  undergraduate  work  in  Electrical  and  Computer  Engi¬ 
neering  in  Taiwan.  William  has  ctmsiderable  parallel  processing  and  super  computing  experience. 
His  dissertation  topic  at  NC  State  University  involved  parallel  Lanczos  algorithms  applied  to 
problems  in  optimization  and  engineering.  He  was  awarded  a  Post  D(x:toral  positicm  with  tiie  U.S. 
Army  High  Performance  Computing  Research  Center  in  Minneapolis,  for  the  academic  years 
1991-92  and  1992-93.  He  accepted  a  position  in  Conq>uter  Science  at  a  Taiwan  university  in 
August  1993. 

•  Sherri  Braxton.  Ms.  Braxton  is  a  very  bright  Afiro-American  U.  S.  student  who  received  her 
undergraduate  degree  with  a  joint  Mathematics  and  Computer  Science  major  at  Wake  Forest 
University  in  May  1992.  Her  participation  in  this  AFOSR  research  project  began  with  her  work  the 
summer  of  1992  on  ill-posed  inverse  problems  in  image  processing.  She  has  a  fellowship  in  the 
PhJ3.  program  in  Computer  Science  at  NC  State  University,  where  the  principle  investigator  is 
directing  her  work. 
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5.  Technical  Publications 


1.  Block  cyclic  SOR  for  Markov  chains  with  p-cyclic  irfinitesinuU  generator,  Lin.  Alg. 

Applications  Special  Issue  on  Iterative  Methods,  vols.  154-156  (1991),  145-224.  (This  79 
page  paper  is  joint  with  K.  Kontovasilis  and  W.  J.  Stewart). 

2.  Adaptive  Lanczos  methods  for  recursive  condition  estimation.  Numerical  Algorithms, 

1  (1991),  1-20.  (with  William  R.  Femg  and  Gene  R  Golub). 

3.  Order-reducing  coryugate  gradients  vs  block  AOR  for  constrained  least  squares  problems,  Lin. 

Alg.  Applications  154-156  (1991),  23-44.  (by  former  Ph.D.  student  AF  Major  Douglas 
James,  and  part  of  his  dissertation  woric  support  by  a  previous  AFOSR  grant). 

4.  Some  fast  Toeplitz  least  squares  algorithms,  Proc.  SPIE  Symposium  on  Signal  Proc. 

Algs.,  Architectures,  and  Implementations,  1566  (1991),  35-46.  (with  J.  Nagy). 

5.  An  inverse  factorization  algorithm  for  linear  prediction,  Lin.  Alg.  Applications  172  (1992), 

169-195  (with  J.  Nagy). 

6.  Fast  adaptive  condition  estimation,  SIAM  J.  Matrix  Anal,  and  Applic.,  13  (1992), 

274-291  (with  D.  Pierce). 

7.  A  parallel  implementation  of  the  inverse  QR  adaptive  filter.  Computers  and  Electrical 

Engineering  18  (1992) ,  291-300  (with  S.  Alexander  and  A.  Ghirinikar). 

8.  Tracking  the  concUthn  number  for  RLS  in  signal  processing.  Mathematics  of  Control, 

Signals,  and  Systems  5  (1992),  23-39  (with  D.  Pierce). 

9.  A  fast  algorithm  for  linear  prediction,  in  Re^nt  Advances  in  Mathematical  Theory  of 

Systems,  Control,  Networks  and  Signal  Processing  11,  Ed.  by  H.  Kimura  and  S. 
Kodama,  MTA  Press  Tokyo  (1992),  15-21  (with  J.  Nagy). 

10.  Block  circulant  preconditioners  for  2-D  deconvolution  problems,  Proc  SPIE  Symposium 

on  Advanced  Signal  Processing  Algorithms,  Architectures,  and  Implement., 
1770  (1992),  60-71.  (with  R.  Chan  and  J.  Nagy). 

11.  Analysis  of  p-cyclic  iterations  for  Markov  chains,  to  appear  in  Linear  Algebra,  Markov 

Chains,  and  Queuing  Models,  Ed.  by  C.  Meyer  and  R.  Plemmons,  Volumes  in 
Mathematics  and  Its  Applications,  Springer  Veilag  48  (1993),  1 1 1-124  (with  A.  Hadjidimos). 

12.  Block  RLS  using  row  Householder  reflections,  Lin.  Alg.  and  Applications,  188  (1993), 

31-62.  (with  A.  Bojanczyk  and  J.  Nagy). 

13.  Fast  recursive  least  squares  using  the  FFT,  IMA  Tech.  Rept  982,  Univ.  Minnesota 

(1992),  Proc  Conf.  on  Mathematics  of  Signal  Processing,  Warwick,  England, 
Oxford  Press,  1993. 

14.  Preconditioned  iterative  regularization  methods  for  ill-posed  problems.  Numerical  Linear 

Algebra  and  Scientific  Computing,  de  Gruyter  Press,  Berlin,  (1993),  141-163. 

(vMth  M.  Hanke  and  J.  Nagy). 
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15.  Iterative  image  restoration  using  FFT-based  preconditioners,  Proc.  30th  Allerton  Conf. 
on  ComnLt  Control  and  Computing,  Allerton,  IL,  1993.  (with  J.  Nagy). 

16.  FFT-based  RLS  in  Signal  Processing,  ICASSP-  93,  Minneapolis,  MN,  IEEE  Press, 

Vol.  in  (1993),  571-574. 

17.  FFT-based  preconditioners  for  Toeplitz-block  least  squares,  SIAM  J.  on  Numerical 
Analysis,  30  (1993),  1740-1768.  (with  R.  Chan  and  J.  Nagy). 

18.  Circulant  preconditioned  Toeplitz  least  squares  iterations,  SIAM  J.  Matrix  Analysis  and 
Applic  15  (1994),  80-97.  (with  R.  Chan  and  J.  Nagy). 

19.  Fast  RLS  adaptive  filtering  by  FFT-based  conjugate  gradient  iterations,  to  appear  in  die 

SIAM  J.  on  Scientific  Computing,  (with  M.  Ng). 

20.  Image  Restoration  using  Fast  Fourier  and  Wavelet  Transforms,  invited  paper  at  the 

International  Symposium  on  Substance  Identification  Technologies,  Conference 
on  Sig^  and  Image  Processing  for  Detection  Systems,  Innsbruck,  Austria,  October  (1993). 
(To  be  published  by  the  Eurc^iean  Optical  Socie^,  1994). 

21.  Fast  LMS-Newton  adaptive  filtering  by  FFT-based  coiyugate  gradient  iterations,  submitted  to 

the  IEEE  Trans.,  on  Signal  Processing,  1994.  (with  M.  Ng). 

22.  Optimizing  closed  loop  adaptive  optics  performance  using  multiple  control  bandwidtiis, 

submitted  to  the  J.  Optical  Soc.  Amer.,  1994.  (with  B.  EUerbroek,  C.  Van  Loan  and  N. 
iKtsianis). 

23.  A  Trace  Maximization  Problem  in  Adaptive  Optics,  in  preparation,  1994.  (with  B.  EUerbroek, 

C.  Van  Loan  and  N.  Pitsianis). 

24.  Fast  restoration  of  atmospherically  blurred  images,  in  pieparatitm,  1994.  (with  J.  Nagy  and 

T.  Torgersen). 
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6.  Invited  Research  Presentations 


1.  Parallel  Algorithms  for  Signal  and  Image  Processing,  Plenary  Talk,  SIAM  Southeast 

Region^  Coftference,  CuUowhee,  NC,  April  1991. 

2.  Non  Negativity  and  Convergent  Iterations,  Invited  Talk,  International  Linear  Algebra 

Society  Conference^  DeKalb,  IL,  May  1991. 

3.  Toeplitz  Least  Squares  Preconditioners,  Invited  Talk,  Southeas  ‘a  Mathematics 

Society  Corference^  Hong  Kong,  June  1991. 

4.  A  Systolic  Array  for  Recursive  Least  Squares  Computations,  ijepartment  of 

Computer  Sci.,  Zhongshan  (Sun  Yat-Sen)  Univ.,  Guangzhou,  P.R.  China,  June  1991. 

5.  An  Inverse  Factorization  Algorithm  for  Linear  Prediction,  Inter.  Confer,  on 

Mathematics  of  Networks  Systems  and  Control,  Kobe,  Japan,  June  1991. 

6.  Extended  SOR  Iterations,  NSF  Conference  on  Multigrid  Methods  for  Partial  Differential 

Equations^  Washington,  DC,  June  1991. 

7.  Fast  Algorithms  for  Toeplitz  Least  Squares  Computations,  SPIE  Symposium  on 

Advanced  Signal  Processing  Algorithrns,  Architec.,  and  Impl.,  San  Diego,  July,  1991. 

8.  Dynamic  Condition  Estimation,  Organizer  of  Mini  symposium,  SIAM  Confer,  on 

Applied  Linear  Algebra,  Minneapolis,  MN,  September  199 1  (co-author  of  2  invited  talks). 

9.  Scalable  Parallel  Algorithms  for  Signal  and  Image  Processing,  Colloquium,  North 

Carolina  Super  Computing  Center,  Research  Triangle  Park,  NC,  November  (1991) 

10.  Matrix  Iterative  Analysis  for  Markov  Chains,  IMA  Workshop  on  Solving  Markov 

Chains,  Minneapolis,  MN,  January  1992  (also,  co-<nganizer  of  the  Workshop). 

11.  FFT-Based  Toeplitz  Least  Squares  Problems,  Invited  Talk,  IMA  Workshop  on 

Sparse  and  Structured  Matrix  Problems,  Minneapolis,  MN,  February  1992. 

12.  Block  Toeplitz  Iterations  on  a  Cray  Y-MP,  Invited  Talk,  Permian  Basin 

Super  Computing  Conference,  Odessa,  TX,  March,  1992. 

13.  Block  Toeplitz  Least  Squares  Iterations,  Colloquium,  NC  State  Univ.,  April  1992. 

14.  Preconditioned  Iterative  Least  Squares  FIR  System  Identification,  Invited  Talk, 

IMA  Workshop  on  Signal  Processing,  Minneapolis,  MN,  April  1992. 

15.  Large-Scale  Block  Toeplitz  Least  Squares  Iterations,  Invited  Talk,  International 

Domaui  Decomposition  for  PDE  Conference,  Conoo,  Italy,  June,  1992. 

16.  Building  FFT-Based  Preconditioners  in  Si^al  Processing,  Invited  Talk, 

Workshop  on  Mathematics  in  Signal  Processing,  Stoors,  CT,  July,  1992. 
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17.  Block  Circulant  Preconditioners  for  2-D  Deconvolution,  SPIE  Symposium  on 

Advanced  Signal  Processing  Algorithms,  Architectures,  and  Implementations,  San  Diego, 
CA,  July,  1^2  (also,  co-organizer  of  Symposium). 

18.  FFT-Based  RLS  in  Signal  Processing,  Invited  Ta^  IMACS  International  Symposium 

on  Mathematical  Moling  and  Scientific*Computing,  Bangalore,  India,  Dec.,  1992. 

19.  Generalized  inverses  in  Signal  Processing,  Invited  Workshop  on  Generalized 

Matrix  Inverses,  Indian  Statistical  Institute,  New  Delhi,  India,  Dec.,  1992. 

20.  Toeplitz  Matrices  and  Iterative  Least  Squares  Estinuition,  Short  Course,  given  at 

the  Indian  Statistical  Institute,  New  Delhi,  India,  Dec.,  1992. 

21.  Fast  Recursive  Least  Squares  using  the  FFT,  Third  IMA  Cotrf.  on  Mathematics  in 

Signal  Processing,  Warwick,  England,  Dec.,  1992. 

22.  FFT-Based  RLS  in  Signal  Processing,  IEEE  Inter.  Cortf.  on  Acoustics,  Speech  and 

Sig.  Processing,  Mimeapolis,  MN,  April,  1993. 

23.  Fast  Recursive  Least  Squares  for  Noise  Cancellation,  Colloquium  Presentation, 

Hiillips  Air  Force  Laboratory,  June,  1993. 

24.  Preconditioners  for  Least  Squares  Iterations,  Plenary  Talk,  Householder  Cortference 

on  Numerical  Analysis,  UCXA  Conf.  Center,  CA,  June  1993. 

25.  Structured  Matrix  Computations  in  Signal  Processing,  Invited  Talk  at  the  SIAM 

Tutorial  on  Numerical  Meth.  in  Control,  Signal  and  Image  Proc.,  Seattle,  August,  1993. 

26.  Building  Preconditioners  for  Toeplitz  and  Related  Systems,  SIAM  Cortf,  on  Lin. 

Alg.  in  Control,  Signal  and  Image  Processing,  Seattle,  WA,  August,  1993. 

27.  Image  Restoration  using  Fast  Fourier  and  Wavelet  Transforms,  Invited  Talk, 

International  Symposium  on  Substance  Identification  Technologies,  Signal  and  Image 
Processing  for  Detection  Systems,  bmsbrucl^  Austria,  October,  1^3. 

28.  An  Optimization  Problem  in  Adaptive  Optics,  Invited  Plenary  Talk,  Lanczos 

Interruttional  Cortference  on  Computational  Mathematics  cmd  Physics,  Raleigh,  NC, 
December,  1993. 
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